PROGRAM MAPS
USE POLYMORPHIC_COMPLEXTAYLOR
TYPE(DAMAP) ROT,SEXT,TOTAL_DA_MAP
REAL(DP) ANGLE,KICK_STRENGTH
type(tree) map_track,map_track_da
type(genfield) genmap
type(pbfield) h
real(dp) ray(lnv),ray_da(lnv),ray2(lnv),ray2_da(lnv),raygen(lnv)
integer mf,me
mf=20
me=21
open(unit=mf,file='results.txt')
open(unit=me,file='results_exact.txt')
CALL INIT(NO1=3,ND1=1,NP1=0,NDPT1 =0)     !   <------------------ init for maps in ND1 degrees of freedom

CALL ALLOC(ROT,SEXT,TOTAL_DA_MAP)
call alloc(map_track,map_track_da); 
call alloc(genmap)
ANGLE=31.0_DP * PI/180.0_DP
KICK_STRENGTH=3.0_DP

call alloc(h)

! This illustrates generating function tracking versus plain Taylor series
! In FPP generating function tracking always neglects the constant part of the map

h=5.01d0*((-angle/2.d0)*((1.d0.mono.'2')+(1.d0.mono.'02')) + KICK_STRENGTH*(1.d0.mono.'3')/3.d0)
TOTAL_DA_MAP=1
TOTAL_DA_MAP=texp(h,TOTAL_DA_MAP)

! Total_map represent a rotation of 5.01*31 degrees modified by a cubic (sextupole-like) term. 

map_track=TOTAL_DA_MAP

!genmap%linear_in=.true.
genmap%ifac=10
genmap=TOTAL_DA_MAP

raygen=0.d0;raygen(1)=0.1d0;
ray=0.d0;ray(1)=0.1d0;

write(mf,'(1x,4(1x,F7.5))') ray(1:2),raygen(1:2)

! Here we track Total_map and the correspounding generating function
do i=1,600
ray=map_track*ray
raygen=genmap*raygen
write(mf,'(1x,4(1x,F7.5))') ray(1:2),raygen(1:2)
enddo

! Here we track one hundreth of map which actually approximate the exact solution in this case.
h=h%h/100.d0
TOTAL_DA_MAP=1
TOTAL_DA_MAP=texp(h,TOTAL_DA_MAP)

ray=0.d0;ray(1)=0.1d0;

map_track=TOTAL_DA_MAP
write(me,'(1x,2(1x,F7.5))') ray(1:2) 
do i=1,60000
ray=map_track*ray
if(mod(i,100)==0) write(me,'(1x,2(1x,F7.5))') ray(1:2)
enddo


call KILL(h)
call KILL(genmap)
CALL KILL(ROT,SEXT,TOTAL_DA_MAP)
call KILL(map_track,map_track_da);
close(mf);close(me)
END PROGRAM MAPS

